Statistical investigation for the article ‘The performance of L2 French children on the LITMUS-QU Nonword repetition task during their first year of exposure: impact of age, vocabulary size, verbal-short term memory and phonological awareness’

---

1 Preparatory steps for the Markdown file

Setting the options for knitr…

library(knitr)
knitr::opts_chunk$set(echo = TRUE,
                       comment = NA,
                       prompt  = FALSE,
                       cache   = FALSE,
                       warning = FALSE,
                       message = FALSE,
                       fig.align="center",
                       fig.width = 8.125,
                       out.width = "100%",
                       fig.path = "D:/figures/eL2/150dpi/",
                       dev = c('png', 'tiff'),
                       dev.args = list(png = list(type = "cairo"), tiff = list(compression = 'lzw')),
                       dpi = 150,
                       cache = TRUE,
                       cache.path = "D:/cache/eL2/150dpi/",
                       autodep = TRUE)
options(width = 1000, scipen = 999, knitr.kable.NA = '')

Setting the options for summarytools…

library(summarytools)
st_options(plain.ascii        = FALSE,        # Always use this option in Rmd documents
            style             = "rmarkdown", # Always use this option in Rmd documents
            footnote          = NA,          # Makes html-rendered results more concise
            subtitle.emphasis = FALSE)  # Improves layout with some rmarkdown themes
st_css()

Loading libraries…

library(tidyverse)
library(ggplot2) # graphics
library(lme4) # glmer
library(lmerTest) # p-values in the summary() and anova() functions
library(performance) # Assess multicolinearity in regression models
library(gridExtra)
library(kableExtra)
library(emmeans) # multiple comparisons and plotting effects
library(car) # For Levene's test
library(FSA) # For Dunn's test
library(stats) # For Fligner-Killeen's test
library(DHARMa) # simulate residuals
library(sjPlot) # Plotting effects
library(glmmTMB)
library(stargazer)

Cleaning the memory…

rm(list = ls(all.names = TRUE))

Specifying a seed for random number generation, for the sake of reproducibility:

set.seed(123)

2 Preparing the data

2.1 Main dataset

Loading the data and building a tibble…

df <- read.table("data.txt", sep = "\t", dec = ".", header = TRUE)
df <- as_tibble(df)

Transforming categorical variables (code, L1, nonword, sex, coda_l, final_l, branching_onset, repetition) into factors…

df <- df %>% mutate(code = as.factor(code),
                    L1 = as.factor(L1),
                    nonword = as.factor(nonword),
                    sex = as.factor(sex),
                    coda_l = as.factor(coda_l),
                    final_l = as.factor(final_l),
                    branching_onset = as.factor(branching_onset),
                    repetition = as.factor(repetition))

Given that we only have 1, 2 or 3 vowels for the nonwords, we will consider V as a categorical variable and therefore recode it as a factor…

df <- df %>% mutate(V = as.factor(V))

We also rename the column code into subject and turn it into a factor…

df <- df %>%
  rename(subject = code) %>%
  mutate(subject = as.factor(subject))

We need to be careful and note that although we have two categorical variables coda_l (whether l appears in internal coda position in the nonword) and final_l (whether l appears in final position in the nonword), we do not have 4 different configurations. Indeed, there is no word in which l appears both in internal coda position and in final position, as shown below:

my_table <- df %>%
  select(nonword, coda_l, final_l) %>%
  unique() %>%
  with(., table(coda_l, final_l)) %>%
  as.data.frame()
colnames(my_table) <- c("l in internal coda position", "l in final position", "# nonwords")
my_table

  l in internal coda position l in final position # nonwords
1                           0                   0         60
2                           1                   0          4
3                           0                   1          7
4                           1                   1          0

Rather than working with coda_l and final_l, we compute a new categorical variable, occurrence_l, with 3 levels: coda for l appearing in internal coda position, final for l appearing in final position and other for nonwords without l in final or coda position (either with l in another position or without l):

df <- df %>% mutate(occurrence_l = if_else(coda_l == 1, "coda", if_else(final_l == 1, "final", "other")),
                    occurrence_l = as.factor(occurrence_l))

We can further observe the distribution of nonwords with respect to the intersections of the variables occurrence_l, V and branching_onset (the number of branching onsets in a nonword).

We first create a data table for nonwords:

df_nonwords <- df %>% select(nonword, occurrence_l, coda_l, final_l, branching_onset, V) %>% unique()

We then inspect several configurations:

df_nonwords %>% with(., table(occurrence_l, branching_onset, dnn = list("occurrence of l", "  branching onset")))

                 branching onset
occurrence of l  0  1  2
          coda   4  0  0
          final  5  2  0
          other 33 25  2

df_nonwords %>% with(., table(occurrence_l, V, dnn = list("occurrence of l", "  # of vowels")))

                 # of vowels
occurrence of l  1  2  3
          coda   0  2  2
          final  3  2  2
          other 23 21 16

df_nonwords %>% with(., table(V, branching_onset, dnn = list("# of vowels", "  branching onset")))

             branching onset
# of vowels  0  1  2
          1 12 14  0
          2 16  7  2
          3 14  6  0

2.2 Additional data about the subjects’s performance

To draw some of the figures of the paper, we need an additional table of information about the subjects’ performances. These data are participant-based rather than item-based.

df_add_data <- read.table("additional_data.txt", header = TRUE, sep = "\t") %>% as_tibble()
summary(df_add_data)

 Participant         NWR_pc_Total      Coda_l_pc        codas_pc      Branching_onset_pc     sC_pc          Final_l_pc          PVC              PCC       
 Length:62          Min.   :0.4366   Min.   :0.000   Min.   :0.5625   Min.   :0.5161     Min.   :0.0000   Min.   :0.0000   Min.   :  9.00   Min.   :82.00  
 Class :character   1st Qu.:0.7077   1st Qu.:0.500   1st Qu.:0.8750   1st Qu.:0.8065     1st Qu.:1.0000   1st Qu.:0.5714   1st Qu.: 93.00   1st Qu.:88.00  
 Mode  :character   Median :0.8169   Median :0.875   Median :0.9375   Median :0.8548     Median :1.0000   Median :0.8571   Median : 95.00   Median :94.00  
                    Mean   :0.7846   Mean   :0.754   Mean   :0.9073   Mean   :0.8460     Mean   :0.9552   Mean   :0.7258   Mean   : 93.23   Mean   :92.49  
                    3rd Qu.:0.8732   3rd Qu.:1.000   3rd Qu.:1.0000   3rd Qu.:0.9355     3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.: 97.00   3rd Qu.:96.00  
                    Max.   :0.9577   Max.   :1.000   Max.   :1.0000   Max.   :1.0000     Max.   :1.0000   Max.   :1.0000   Max.   :100.00   Max.   :99.00  
                                                                                                                                            NA's   :1      

(pc means percentage correct)

2.3 Additional data about L1 syllabic complexity

The data introduced in this subsection come from the Lapsyd database (Maddieson et al., 2013, Proc. of the 14th Annual Conf. of the Intl. Speech Communication Association (Interspeech 2013), pp. 3022-3026). Among other features, this database describes, for most languages, the complexity of the onset, the nucleus and the coda of their syllables with scores from 0 (least complex) to 3 (most complex), with 1 being the default score.

Loading the data and building a tibble…

df_L1 <- read.table("Data on the syllabic complexity of L1.txt", sep = "\t", dec = ".", header = TRUE, stringsAsFactors = T) %>% as_tibble()

We have data regarding the complexity of the syllables for all the languages of our study but Kabyle:

df_L1

# A tibble: 17 × 2
   Language             L1_syll_complexity
   <fct>                             <int>
 1 Arabic                                4
 2 Arabic/Italian                        6
 3 Italian                               6
 4 Portuguese                            5
 5 Mandarin                              4
 6 Ukrainian/Russian                     7
 7 Romanian                              7
 8 Albanian                              7
 9 Brazilian Portuguese                  5
10 Armenian                              6
11 Russian/Armenian                      7
12 Japanese                              4
13 Georgian/German                       8
14 Russian/Chechen                       8
15 Russian                               7
16 Arabic/Spanish                        5
17 Spanish                               5

We integrate these data to our main dataset:

df <- df %>% left_join(df_L1, by = c("L1" = "Language"))

3 Description of the problems and of the data

3.1 Description of the investigation and of the variables

Our experiment consists in presenting French nonwords to children acquiring this language as an L2, and assess their ability to repeat these target nonwords. We investigate which parameters may explain success or failure at this task. To this end, our statistical approach rests on a (valid) mixed-effects logistic regression model.

The dependent variable of our model is repetition (0 for failure to repeat, 1 for success)

As for the meaningful independent variables, we consider:

• occurrence_l (whether l appears in coda or final position) (‘coda’ if l appears in internal coda position, ‘final’ if l appears in final position, ‘other’ for other configurations)
• branching_onset (number of branching onsets in the nonword) (0, 1 or 2)
• V (number of vowels of the nonword) (1, 2 or 3)
• rec_voc (size of a child’s French receptive vocabulary)
• phono_mem (size of a child’s phonological memory)
• phono_awareness (z-scored score reflecting the child’s phonological awareness)
• L1_syll_complexity (overall complexity of the syllabic structures in the child’s L1)
• LoE (length of a child’s exposure to French, in days)
• age (child’s age, in months)

We are not only interested in main (fixed) effects, but also in the possible interactions between the previous predictors. We do not consider triple (or even higher-order) interactions for the sake of clarity, and thus focus on simple (i.e., ‘double’ interactions). Even with this approach, however, there are 21 possible interactions given that we have 7 primary predictors (7!/(5!2!)). Rather than considering all possible options and get lost in backward and forward stepwise regression, we choose to investigate only a number of ‘target’ interactions, which we deem meaningful for our study and given our expertise in the domain. These interactions are the following:

• occurrence_l * LoE
• branching_onset * LoE
• rec_voc * LoE
• age * phono_mem
• age * rec_voc
• phono_mem * rec_voc

We finally consider a number of random effects to capture the non-independence of our observations:

• subject (with 62 different levels - one per subject)
• L1 (with 18 different levels, covering both monolingual and bilingual subjects)
• nonword (with 71 different levels - one per nonword)

subject is nested into L1, but since there is no ambiguity (i.e. two subjects with the same name/code speaking different languages), nesting subject into L1 explicitly in the model, i.e., (1|L1/subject), is similar to considering crossed random effects, i.e. (1|L1) + (1|subject).

We can additionally formulate precise oriented hypotheses to be tested for the main effects (all other things being equal):

• Regarding occurrence_l: a nonword gets easier to repeat when shifting from l in coda position to l in final position to another structure
• Regarding branching_onset: the more branching onsets a nonword has, the more difficult it is to repeat
• Regarding V: following the litterature, a nonword becomes more difficult to repeat when it has 3 vowels (no difference between 1 and 2 vowels)
• Regarding, rec_voc: the larger a child’s French receptive vocabulary, the easier it is for them to repeat the nonwords
• Regarding phono_mem: the larger a child’s phonological memory, the easier it is for them to repeat the nonwords
• Regarding phono_awareness: the more developed the children’s phonological awareness, the easier it is for them to repeat the nonwords
• Regarding L1_syll_complexity: the more complex the syllables can be in the children’s L1, the easier it is for them to repeat the nonwords in French
• Regarding LoE: the longer the exposure to French, the easier it is for them to repeat the nonwords
• Regarding age: the older the children, the easier it is for them to repeat the nonwords

The hypotheses for Onset_complexity and Coda_complexity derive from the fact that both onset and coda can be complex in French - with scores of 3 for both of them in Lapsyd.

3.2 Overview of the observations

There are 4402 observations.

We can inspect the distributions of values for the different variables:

df %>%
  dfSummary() %>%
  print(method = 'render', footnote = NA) 

Data Frame Summary

df
Dimensions: 4402 x 18
Duplicates: 0

No	Variable	Stats / Values	Freqs (% of Valid)	Graph	Valid	Missing
1	subject [factor]	1. EANA-A1 2. EANA-A2 3. EANA-A3 4. EANA-A4 5. EANA-A5 6. EANA-A6 7. EANA-A7 8. EANA-A8 9. EANA-A9 10. EANA-AL1 [ 52 others ]	71 ( 1.6% ) 71 ( 1.6% ) 71 ( 1.6% ) 71 ( 1.6% ) 71 ( 1.6% ) 71 ( 1.6% ) 71 ( 1.6% ) 71 ( 1.6% ) 71 ( 1.6% ) 71 ( 1.6% ) 3692 ( 83.9% )		4402 (100.0%)	0 (0.0%)
2	sex [factor]	1. F 2. M	1704 ( 38.7% ) 2698 ( 61.3% )		4402 (100.0%)	0 (0.0%)
3	L1 [factor]	1. Albanian 2. Arabic 3. Arabic/Italian 4. Arabic/Spanish 5. Armenian 6. Brazilian Portuguese 7. Georgian/German 8. Italian 9. Japanese 10. Kabyle [ 8 others ]	426 ( 9.7% ) 639 ( 14.5% ) 852 ( 19.4% ) 71 ( 1.6% ) 213 ( 4.8% ) 284 ( 6.5% ) 71 ( 1.6% ) 355 ( 8.1% ) 142 ( 3.2% ) 71 ( 1.6% ) 1278 ( 29.0% )		4402 (100.0%)	0 (0.0%)
4	age [integer]	Mean (sd) : 89.6 (10) min ≤ med ≤ max: 72 ≤ 89.5 ≤ 109 IQR (CV) : 17 (0.1)	31 distinct values		4402 (100.0%)	0 (0.0%)
5	LoE [integer]	Mean (sd) : 180 (79.3) min ≤ med ≤ max: 24 ≤ 180 ≤ 350 IQR (CV) : 118 (0.4)	44 distinct values		4402 (100.0%)	0 (0.0%)
6	nonword [factor]	1. faku 2. fal 3. fapul 4. fapus 5. fikapul 6. fikupla 7. fikuspa 8. filpa 9. fips 10. fiska [ 61 others ]	62 ( 1.4% ) 62 ( 1.4% ) 62 ( 1.4% ) 62 ( 1.4% ) 62 ( 1.4% ) 62 ( 1.4% ) 62 ( 1.4% ) 62 ( 1.4% ) 62 ( 1.4% ) 62 ( 1.4% ) 3782 ( 85.9% )		4402 (100.0%)	0 (0.0%)
7	V [factor]	1. 1 2. 2 3. 3	1612 ( 36.6% ) 1550 ( 35.2% ) 1240 ( 28.2% )		4402 (100.0%)	0 (0.0%)
8	coda_l [factor]	1. 0 2. 1	4154 ( 94.4% ) 248 ( 5.6% )		4402 (100.0%)	0 (0.0%)
9	branching_onset [factor]	1. 0 2. 1 3. 2	2604 ( 59.2% ) 1674 ( 38.0% ) 124 ( 2.8% )		4402 (100.0%)	0 (0.0%)
10	final_l [factor]	1. 0 2. 1	3968 ( 90.1% ) 434 ( 9.9% )		4402 (100.0%)	0 (0.0%)
11	repetition [factor]	1. 0 2. 1	948 ( 21.5% ) 3454 ( 78.5% )		4402 (100.0%)	0 (0.0%)
12	rec_voc [integer]	Mean (sd) : 19.7 (6.4) min ≤ med ≤ max: 5 ≤ 20 ≤ 31 IQR (CV) : 10 (0.3)	23 distinct values		4402 (100.0%)	0 (0.0%)
13	phono_mem [integer]	Mean (sd) : 3.8 (0.7) min ≤ med ≤ max: 2 ≤ 4 ≤ 6 IQR (CV) : 1 (0.2)	2 : 142 ( 3.2% ) 3 : 1136 ( 25.8% ) 4 : 2485 ( 56.5% ) 5 : 568 ( 12.9% ) 6 : 71 ( 1.6% )		4402 (100.0%)	0 (0.0%)
14	rec_voc_zscore [numeric]	Mean (sd) : -13.8 (6.1) min ≤ med ≤ max: -30.8 ≤ -13.4 ≤ -3.6 IQR (CV) : 8.3 (-0.4)	38 distinct values		4402 (100.0%)	0 (0.0%)
15	phono_mem_zscore [numeric]	Mean (sd) : -1.2 (0.9) min ≤ med ≤ max: -3.3 ≤ -1 ≤ 1.5 IQR (CV) : 1.2 (-0.8)	11 distinct values		4402 (100.0%)	0 (0.0%)
16	phono_awareness [numeric]	Mean (sd) : -2.8 (1.5) min ≤ med ≤ max: -4 ≤ -3.9 ≤ 0.6 IQR (CV) : 1.7 (-0.5)	17 distinct values		4047 (91.9%)	355 (8.1%)
17	occurrence_l [factor]	1. coda 2. final 3. other	248 ( 5.6% ) 434 ( 9.9% ) 3720 ( 84.5% )		4402 (100.0%)	0 (0.0%)
18	L1_syll_complexity [integer]	Mean (sd) : 5.7 (1.1) min ≤ med ≤ max: 4 ≤ 6 ≤ 8 IQR (CV) : 2 (0.2)	4 : 852 ( 19.7% ) 5 : 923 ( 21.3% ) 6 : 1420 ( 32.8% ) 7 : 994 ( 23.0% ) 8 : 142 ( 3.3% )		4331 (98.4%)	71 (1.6%)

4 First figures and statistical analyses

4.1 Boxplot for children’s percentage of correct repetition

df_add_data %>%
  mutate(group = "") %>%
  ggplot(aes(x = group, y = NWR_pc_Total)) + 
  geom_boxplot(outlier.size = 3, fill="#A4A4A4") +
  stat_summary(fun = mean, geom = "point", shape = 4, size = 4, stroke = 2) +
  geom_jitter(shape = 16, size = 1, position = position_jitter(0.075)) +
  theme_classic() +
  ylab("Percentage of correct repetition") + xlab("") +
  scale_y_continuous(labels = scales::percent, breaks = seq(0.0, 1.0, 0.05))

4.2 Boxplot for receptive vocabulary size (z score) (by subject)

df %>%
  group_by(subject) %>%
  summarize(rec_voc_zscore = mean(rec_voc_zscore)) %>%
  ungroup() %>%
  mutate(group = as.factor("participants")) %>%
  ggplot(aes(x = group, y = rec_voc_zscore)) + 
  geom_boxplot(outlier.size = 3, fill="#A4A4A4", coef = 1.5) +
  stat_summary(fun = mean, geom = "point", shape = 4, size = 4, stroke = 2) +
  geom_jitter(shape = 16, size = 1, position = position_jitter(0.075)) +
  theme_classic() +
  ylab("Receptive vocabulary size (z score)") + xlab("") +
  scale_y_continuous(breaks = seq(-30.0, 0.0, 5.0))

4.3 Boxplot for the percentage of correct repetition (by subject and by 5 types of syllable structures)

We first prepare the data for the figure…

df_struct <- df_add_data %>%
  select(-PCC, -PVC, -NWR_pc_Total) %>%
  rename(`Coda /l/` = Coda_l_pc, `Branching onset` = Branching_onset_pc, `Final /l/` = Final_l_pc, `Obstruent coda` = codas_pc, `#sC cluster` = sC_pc) %>%
  pivot_longer(cols = 2:6, names_to = "struct", values_to = "pc_correct_rep") %>%
  mutate(struct = as.factor(struct),
         struct = fct_relevel(struct, "Branching onset", "Obstruent coda", "Coda /l/", "Final /l/", "#sC cluster"))

Then we draw the boxplot…

df_struct %>%
  ggplot(aes(x = struct, y = pc_correct_rep, fill = struct)) + 
  geom_boxplot(outlier.size = 3, show.legend = FALSE, colour = "black") +
  scale_fill_grey() +
  stat_summary(fun = mean, geom = "point", shape = 4, size = 4, stroke = 2, show.legend = FALSE) +
  geom_jitter(shape = 16, size = 1, position = position_jitter(0.075), show.legend = FALSE) +
  theme_classic() +
  ylab("Percentage of correct repetition") + xlab("") +
  scale_y_continuous(labels = scales::percent, breaks = seq(0.0, 1.0, 0.1))

4.4 Comparison between the median values of the percentages of correct repetition for the 5 types of syllable structures

In this section and the following one, given the non-normality of the distribution of percentages of correct repetition in our 5 groups, we need to rely on non-parametric tests.

We apply a Kruskal-Wallis test to check whether there are differences in terms of median between the 5 groups:

df_struct %>% kruskal.test(pc_correct_rep ~ struct, data = .)

    Kruskal-Wallis rank sum test

data:  pc_correct_rep by struct
Kruskal-Wallis chi-squared = 59.564, df = 4, p-value = 0.000000000003582

We can clearly reject the null hypothesis of no difference between the medians.

df_struct %>% group_by(struct) %>% summarize(median = median(pc_correct_rep)) %>% arrange(-median)

# A tibble: 5 × 2
  struct          median
  <fct>            <dbl>
1 #sC cluster      1    
2 Obstruent coda   0.938
3 Coda /l/         0.875
4 Final /l/        0.857
5 Branching onset  0.855

We run post-hoc tests to look at pair differences:

df_struct %>% dunnTest(pc_correct_rep ~ struct, data = ., method = "holm")

                         Comparison           Z            P.unadj             P.adj
1     #sC cluster - Branching onset  6.47529641 0.0000000000946260 0.000000000946260
2            #sC cluster - Coda /l/  5.06470605 0.0000004090309978 0.000003272247983
3        Branching onset - Coda /l/ -1.41059035 0.1583654377394939 0.475096313218482
4           #sC cluster - Final /l/  6.39404724 0.0000000001615515 0.000000001453963
5       Branching onset - Final /l/ -0.08124917 0.9352438002034137 0.935243800203414
6              Coda /l/ - Final /l/  1.32934119 0.1837354315793458 0.367470863158692
7      #sC cluster - Obstruent coda  3.01513033 0.0025686885456720 0.012843442728360
8  Branching onset - Obstruent coda -3.46016608 0.0005398423202431 0.003778896241702
9         Coda /l/ - Obstruent coda -2.04957573 0.0404058511140099 0.161623404456040
10       Final /l/ - Obstruent coda -3.37891691 0.0007277199779395 0.004366319867637

We observe clear differences between (the median of % of correct repetitions for) “#sC cluster” and the 4 other groups. We also find a significant difference between “Branching onset” and “Obstruent coda”, as well as between “Final /l/” and “Obstruent Coda’.

4.5 Comparison between the variances of the percentages of correct repetition for the 5 types of syllable structures

We can assess the homogeneity of the variance for the five types of syllable structures with Fligner-Killeen (median) test (which is non-parametric, unlike Levene’s test):

df_struct %>% fligner.test(pc_correct_rep ~ struct, data = .)

    Fligner-Killeen test of homogeneity of variances

data:  pc_correct_rep by struct
Fligner-Killeen:med chi-squared = 108.58, df = 4, p-value < 0.00000000000000022

The test clearly reject the null hypothesis of equal variances in the 5 groups.

df_struct %>% group_by(struct) %>% summarize(var = var(pc_correct_rep)) %>% arrange(-var)

# A tibble: 5 × 2
  struct             var
  <fct>            <dbl>
1 Coda /l/        0.0953
2 Final /l/       0.0818
3 #sC cluster     0.0233
4 Obstruent coda  0.0137
5 Branching onset 0.0122

The 5 variances may be arranged into two groups: the variance is quite clearly higher for “Coda /l/” and “Final /l/” than for “#sC cluster”, “Obstruent coda” and “Branching onset.

We can look at pair differences between the variances of the 5 groups:

pairs <- c(
  "Coda /l/ - Final /l/",
  "Coda /l/ - #sC cluster",
  "Coda /l/ - Obstruent coda",
  "Coda /l/ - Branching onset",
  "Final /l/ - #sC cluster",
  "Final /l/ - Obstruent coda",
  "Final /l/ - Branching onset",
  "#sC cluster - Obstruent coda",
  "#sC cluster - Branching onset",
  "Obstruent coda - Branching onset"
)

test_statistics <- c(
  df_struct %>% filter(struct %in% c("Coda /l/", "Final /l/")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("Coda /l/", "#sC cluster")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("Coda /l/", "Obstruent coda")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("Coda /l/", "Branching onset")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("Final /l/", "#sC cluster")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("Final /l/", "Obstruent coda")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("Final /l/", "Branching onset")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("#sC cluster", "Obstruent coda")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("#sC cluster", "Branching onset")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic,
  df_struct %>% filter(struct %in% c("Obstruent coda", "Branching onset")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$statistic
)

p_values <- c(
  df_struct %>% filter(struct %in% c("Coda /l/", "Final /l/")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("Coda /l/", "#sC cluster")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("Coda /l/", "Obstruent coda")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("Coda /l/", "Branching onset")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("Final /l/", "#sC cluster")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("Final /l/", "Obstruent coda")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("Final /l/", "Branching onset")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("#sC cluster", "Obstruent coda")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("#sC cluster", "Branching onset")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value,
  df_struct %>% filter(struct %in% c("Obstruent coda", "Branching onset")) %>% fligner.test(pc_correct_rep ~ struct, data = .) %>% .$p.value
)

variance_pair_differences <- tibble(pair = pairs, statistic = test_statistics, p_value = p_values) %>%
  mutate(adjusted_p_values = p.adjust(p_values, method = "holm"))

variance_pair_differences

# A tibble: 10 × 4
   pair                             statistic  p_value adjusted_p_values
   <chr>                                <dbl>    <dbl>             <dbl>
 1 Coda /l/ - Final /l/                 3.46  6.28e- 2          1.26e- 1
 2 Coda /l/ - #sC cluster              53.7   2.39e-13          2.39e-12
 3 Coda /l/ - Obstruent coda           45.8   1.29e-11          1.16e-10
 4 Coda /l/ - Branching onset          45.4   1.59e-11          1.28e-10
 5 Final /l/ - #sC cluster             39.9   2.63e-10          1.84e- 9
 6 Final /l/ - Obstruent coda          30.5   3.35e- 8          1.67e- 7
 7 Final /l/ - Branching onset         25.0   5.81e- 7          2.32e- 6
 8 #sC cluster - Obstruent coda        24.1   9.26e- 7          2.78e- 6
 9 #sC cluster - Branching onset       32.5   1.18e- 8          7.06e- 8
10 Obstruent coda - Branching onset     0.105 7.46e- 1          7.46e- 1

All differences are significant except between “Coda /l/” and “Final /l/”, and between “Obstruent coda” and “Branching onset”.

5 Building a valid regression model

5.1 Preparing a ‘sub-dataset’ for the regression model

We have a few missing data for L1_syll_complexity and phono_awareness. Given these are variables of interest and we want to include them as predictors in our model, we need to extract a subset of observations.

We first drop a participant speaking Kabyle, given that we don’t have data about the syllabic complexity in this language:

df_reduced <- df %>% filter(! is.na(L1_syll_complexity))

Second, a few participants did not understand the task assessing their phonological awareness and therefore do not have any related score:

df_reduced %>% filter(is.na(phono_awareness)) %>% group_by(L1, subject) %>% tally()

# A tibble: 5 × 3
# Groups:   L1 [3]
  L1       subject      n
  <fct>    <fct>    <int>
1 Albanian EANA-AL4    71
2 Japanese EANA-J1     71
3 Japanese EANA-J2     71
4 Romanian EANA-R1     71
5 Romanian EANA-R2     71

They consist in one Albanian, two Japanese and two Romanian participants. We drop these participants from our final dataset:

df_reduced <- df_reduced %>% filter(! is.na(phono_awareness))

5.2 Defining the parameters for the glmer model.

We first define a control structure for our glmer models…

ctrl <- glmerControl(optimizer="bobyqa", optCtrl = list(maxfun=100000), calc.derivs = T, check.conv.grad = .makeCC("warning", tol = 1e-3, relTol = NULL))

5.3 Computing the model

We define the model with the function glmer()

model_rep <- glmer(repetition ~ occurrence_l + branching_onset + V + scale(LoE) + scale(age) + scale(rec_voc) + scale(phono_mem) +
                       scale(phono_awareness) + L1_syll_complexity +
                       occurrence_l:scale(LoE) +
                       branching_onset:scale(LoE) +
                       scale(rec_voc):scale(LoE) +
                       scale(age):scale(phono_mem) +
                       scale(age):scale(rec_voc) +
                       scale(rec_voc):scale(phono_mem) +
                       (1 | subject) +
                       (1 | L1) + 
                       (1 | nonword),
                     data = df_reduced, control = ctrl, family = binomial(link = "logit"))

Our model converges without any problem. We can confirm it is not singular (wich could affect its poxer, i.e., its capacity to detect true effects (Matuschek et al., 2017)).

Checking the singularity of the model…

isSingular(model_rep)

[1] FALSE

We can look at the summary of the model to inspect whether there is any obvious problem:

summary(model_rep)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: repetition ~ occurrence_l + branching_onset + V + scale(LoE) +      scale(age) + scale(rec_voc) + scale(phono_mem) + scale(phono_awareness) +      L1_syll_complexity + occurrence_l:scale(LoE) + branching_onset:scale(LoE) +      scale(rec_voc):scale(LoE) + scale(age):scale(phono_mem) +      scale(age):scale(rec_voc) + scale(rec_voc):scale(phono_mem) +      (1 | subject) + (1 | L1) + (1 | nonword)
   Data: df_reduced
Control: ctrl

     AIC      BIC   logLik deviance df.resid 
  3616.1   3767.0  -1784.0   3568.1     3952 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.1381  0.2091  0.3445  0.4950  1.7302 

Random effects:
 Groups  Name        Variance Std.Dev.
 nonword (Intercept) 0.2953   0.5434  
 subject (Intercept) 0.2380   0.4878  
 L1      (Intercept) 0.0351   0.1874  
Number of obs: 3976, groups:  nonword, 71; subject, 56; L1, 16

Fixed effects:
                                 Estimate Std. Error z value   Pr(>|z|)    
(Intercept)                      0.804346   0.770716   1.044    0.29665    
occurrence_lfinal                0.351099   0.400355   0.877    0.38050    
occurrence_lother                1.364402   0.336003   4.061 0.00004893 ***
branching_onset1                -0.812214   0.171910  -4.725 0.00000231 ***
branching_onset2                -1.830245   0.464614  -3.939 0.00008173 ***
V2                              -0.101115   0.198948  -0.508    0.61128    
V3                              -0.595412   0.202572  -2.939    0.00329 ** 
scale(LoE)                       0.182079   0.173767   1.048    0.29471    
scale(age)                       0.155004   0.108546   1.428    0.15329    
scale(rec_voc)                   0.158190   0.119882   1.320    0.18699    
scale(phono_mem)                 0.038202   0.092159   0.415    0.67849    
scale(phono_awareness)           0.250440   0.105636   2.371    0.01775 *  
L1_syll_complexity               0.038894   0.115452   0.337    0.73621    
occurrence_lfinal:scale(LoE)     0.135684   0.197825   0.686    0.49279    
occurrence_lother:scale(LoE)    -0.078190   0.167454  -0.467    0.64055    
branching_onset1:scale(LoE)     -0.151053   0.093333  -1.618    0.10557    
branching_onset2:scale(LoE)     -0.009488   0.217244  -0.044    0.96517    
scale(LoE):scale(rec_voc)        0.048142   0.100188   0.481    0.63086    
scale(age):scale(phono_mem)      0.020501   0.089731   0.228    0.81928    
scale(age):scale(rec_voc)       -0.087775   0.081682  -1.075    0.28255    
scale(rec_voc):scale(phono_mem)  0.070399   0.140309   0.502    0.61585    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Everything appears ok at first sight, but we need to further assess the validity of your model.

5.4 Assessing the validity of the model

In order to trust the outputs of our model, we need to make sure that the assumptions underlying a mixed-effects logistic regression are met. While these assumptions are not as clear as for linear regression, the following conditions matter: - The normality of the different distributions for the random effects - The correct distribution of the residuals of the model - The absence of problematic multicollinearity

(Note that we already know that our model is not singular)

5.4.1 Normality of the levels of the random effects

p <- plot_model(model_rep, type = "diag", show.values = TRUE, value.offset = .3)
p$nonword + labs(title = "Random intercepts for nonword")

p$L1 + labs(title = "Random intercepts for L1")

p$subject + labs(title = "Random intercepts for subject") + theme_sjplot()

The quantile-quantile plots of the various random effects indicate that their levels follow a distribution which is quite close to a Gaussian distribution.

5.4.2 Residuals of the model

There is no expected normality of the model residuals with a logistic regression, but we can simulate residuals with the function simulateResiduals() of the DHARMa package. According to the authors of the package: ‘Even experienced statistical analysts currently have few options to diagnose misspecification problems in GLM. DHARMa package aims at solving these problems by creating readily interpretable residuals for generalized linear models that are standardized to values between 0 and 1, and that can be interpreted as intuitively as residuals for the linear model. This is achieved by a simulation-based approach that transforms the residuals to a standardized scale. The key idea is that, if the model is correctly specified, then the observed data should look like as if it was created from the fitted model. Hence, for a correctly specified model, all values of the cumulative distribution should appear with equal probability. That means we expect the distribution of the residuals to be flat, regardless of the binomial model structure.’ — Source: DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models

sim.residuals <- simulateResiduals(fittedModel = model_rep, n = 1000, refit = FALSE)
plotSimulatedResiduals(simulationOutput = sim.residuals)

According to the graphics, the simulated residuals do not reveal any misspecification of the model - they are uniformly and homogeneously distributed with respect to the model predictions (right graphic), and they also appear to be quite linearly distributed.

We can also inspect the residuals not with respect to the model predictions, but to the different main fixed effects.

plotResiduals(sim.residuals, form = df_reduced$occurrence_l, xlab = "Occurrence of l", ylab = "Simulated Standardized Residuals")

plotResiduals(sim.residuals, form = df_reduced$branching_onset, xlab = "branching_onset", ylab = "Simulated Standardized Residuals")

plotResiduals(sim.residuals, form = df_reduced$V, xlab = "V", ylab = "Simulated Standardized Residuals", asFactor = F)

We do not see any violation of the uniformity and homogeneity of the distribution of the simulated residuals for our different categorical predictors.

plotResiduals(sim.residuals, form = df_reduced$rec_voc, xlab = "rec_voc", ylab = "Simulated Standardized Residuals")

plotResiduals(sim.residuals, form = df_reduced$phono_mem, xlab = "phono_mem", ylab = "Simulated Standardized Residuals", asFactor = F)

plotResiduals(sim.residuals, form = df_reduced$LoE, xlab = "LoE", ylab = "Simulated Standardized Residuals")

plotResiduals(sim.residuals, form = df_reduced$age, xlab = "age", ylab = "Simulated Standardized Residuals")

plotResiduals(sim.residuals, form = df_reduced$L1_syll_complexity, asFactor = F, xlab = "Syllabic complexity of the L1", ylab = "Simulated Standardized Residuals")

plotResiduals(sim.residuals, form = df_reduced$phono_awareness, xlab = "Phonological awareness", ylab = "Simulated Standardized Residuals")

We also do not see any violation of the uniformity and homogeneity of the distribution of the simulated residuals for our different numerical predictors.

5.4.3 Multicollinearity

We finally need to check whether there is multicolinearity in the model, i.e., whether the predictive effect of a fixed predictor is not too close to the predictive effect of a linear combination of other fixed predictors. We use the function check_collinearity from the performance package, as it can be applied to mixed-effects logistic regressions:

check_collinearity(model_rep)

# Check for Multicollinearity

Low Correlation

                            Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
                    occurrence_l 1.11 [1.08, 1.15]         1.05      0.90     [0.87, 0.93]
                 branching_onset 1.22 [1.18, 1.26]         1.10      0.82     [0.79, 0.85]
                               V 1.18 [1.14, 1.22]         1.09      0.85     [0.82, 0.88]
                      scale(LoE) 4.68 [4.43, 4.95]         2.16      0.21     [0.20, 0.23]
                      scale(age) 1.78 [1.70, 1.86]         1.33      0.56     [0.54, 0.59]
                  scale(rec_voc) 1.84 [1.77, 1.93]         1.36      0.54     [0.52, 0.57]
                scale(phono_mem) 1.27 [1.22, 1.32]         1.12      0.79     [0.76, 0.82]
          scale(phono_awareness) 1.55 [1.49, 1.62]         1.24      0.65     [0.62, 0.67]
              L1_syll_complexity 1.74 [1.66, 1.82]         1.32      0.58     [0.55, 0.60]
         occurrence_l:scale(LoE) 4.78 [4.53, 5.06]         2.19      0.21     [0.20, 0.22]
      branching_onset:scale(LoE) 1.44 [1.38, 1.50]         1.20      0.70     [0.67, 0.72]
       scale(LoE):scale(rec_voc) 1.42 [1.37, 1.48]         1.19      0.70     [0.68, 0.73]
     scale(age):scale(phono_mem) 1.17 [1.14, 1.22]         1.08      0.85     [0.82, 0.88]
       scale(age):scale(rec_voc) 1.21 [1.17, 1.26]         1.10      0.82     [0.79, 0.85]
 scale(rec_voc):scale(phono_mem) 1.70 [1.63, 1.78]         1.31      0.59     [0.56, 0.61]

None of the fixed effects has a VIF higher than 5. Two terms, LoE and occurrence_l:LoE have a VIF close to 5, but it is actually not unexpected that multicollinearity appears with interaction terms in a model, and it is usually not considered a problem, since the estimate of the interaction term will not be biased. Overall, we do not have any obvious problem of multicollinearity.

5.5 Conclusion about the validity of the model

We did not find any issue with our model. It therefore appears to be valid and we can now asses our hypotheses with it.

6 Assessing our hypotheses

6.1 Analyzing the fixed effects

In the next paragraphs, one should note that the graphical representations report probabilities of correct repetition of a nonword, i.e., in the ‘numerical space’ of the predicted variable. However, the statistical tests of the various effects consider the non-transformed estimates and confidence intervals, i.e. in the numerical space of the linear combination of predictors / independent variables.

We rely on the emmeans package and the computation of estimated marginal means of either levels of factors or linear trends. To adjust the p-values in the case of multiple tests, we rely on a multivariate t distribution (adjust = “mvt”) with the same covariance structure as the estimates, as it is the closest to an exact adjustment (see https://cran.r-project.org/web/packages/emmeans/vignettes/confidence-intervals.html).

6.1.1 Investigating the interactions

We can first investigate the higher-level predictors of the model, that is the different interactions.

• occurrence_l * LoE
• branching_onset * LoE
• rec_voc * LoE
• age * phono_mem
• age * rec_voc
• phono_mem * rec_voc

The first two interactions are between one categorical variable and the continuous variable LoE (occurrence_l : LoE and branching_onset : LoE)

We first define a range of variation for the values of LoE

list.LoE <- list(LoE = seq(min(df_reduced$LoE), max(df_reduced$LoE), by = 1))

For occurrence_l : LoE:

plot_model(model_rep, type = "emm", terms = c("LoE [all]", "occurrence_l"))

emtrends(model_rep, pairwise ~ occurrence_l | LoE, var= "LoE", adjust = "mvt", infer = c(T,T))

$emtrends
LoE = 184:
 occurrence_l LoE.trend      SE  df asymp.LCL asymp.UCL z.ratio p.value
 coda          0.001585 0.00238 Inf -0.003089   0.00626   0.665  0.5063
 final         0.003258 0.00204 Inf -0.000744   0.00726   1.595  0.1106
 other         0.000621 0.00141 Inf -0.002143   0.00339   0.440  0.6597

Results are averaged over the levels of: branching_onset, V 
Confidence level used: 0.95 

$contrasts
LoE = 184:
 contrast       estimate      SE  df asymp.LCL asymp.UCL z.ratio p.value
 coda - final  -0.001673 0.00244 Inf  -0.00735   0.00401  -0.686  0.7663
 coda - other   0.000964 0.00206 Inf  -0.00384   0.00577   0.467  0.8838
 final - other  0.002637 0.00165 Inf  -0.00120   0.00647   1.600  0.2388

Results are averaged over the levels of: branching_onset, V 
Confidence level used: 0.95 
Conf-level adjustment: mvt method for 3 estimates 
P value adjustment: mvt method for 3 tests 

The figures show that the slopes for the levels other and coda are quite close, and are visually quite different from the slope for the third level final. The statistical tests, however, fail to detect a significant difference. The reason for that is possibly the large standard errors for coda and final, which are due in turn to the small number of nonwords which relate to these two levels (4 and 7 respectively, to be compared to 60 nonwords without l in coda or final position):

my_table <- df_reduced %>%
  select(nonword, occurrence_l) %>%
  unique() %>%
  with(., table(occurrence_l)) %>%
  as.data.frame()
colnames(my_table) <- c("Occurrence of l", "# nonwords")
my_table

  Occurrence of l # nonwords
1            coda          4
2           final          7
3           other         60

For branching_onset : LoE:

plot_model(model_rep, type = "emm", terms = c("LoE [all]", "branching_onset"))

emtrends(model_rep, pairwise ~ branching_onset | LoE, var= "LoE", adjust = "mvt", infer = c(T,T))

$emtrends
LoE = 184:
 branching_onset LoE.trend      SE  df asymp.LCL asymp.UCL z.ratio p.value
 0                0.002481 0.00138 Inf -0.000227   0.00519   1.796  0.0725
 1                0.000619 0.00156 Inf -0.002440   0.00368   0.397  0.6917
 2                0.002364 0.00290 Inf -0.003311   0.00804   0.817  0.4142

Results are averaged over the levels of: occurrence_l, V 
Confidence level used: 0.95 

$contrasts
LoE = 184:
 contrast                             estimate      SE  df asymp.LCL asymp.UCL z.ratio p.value
 branching_onset0 - branching_onset1  0.001863 0.00115 Inf  -0.00079   0.00452   1.618  0.2235
 branching_onset0 - branching_onset2  0.000117 0.00268 Inf  -0.00606   0.00629   0.044  0.9989
 branching_onset1 - branching_onset2 -0.001746 0.00266 Inf  -0.00788   0.00439  -0.656  0.7789

Results are averaged over the levels of: occurrence_l, V 
Confidence level used: 0.95 
Conf-level adjustment: mvt method for 3 estimates 
P value adjustment: mvt method for 3 tests 

While the figures show that the slopes for LoE differ according to the value of branching_onset, the statistical tests reveal that these slopes are actually not significantly different from each other.

We then have 4 interactions which consist of two continuous variables : rec_voc : LoE, age : phono_mem, months : rec_voc and phono_mem : rec_voc

For rec_voc : LoE:

plot_model(model_rep, type = "emm", terms = c("LoE [all]", "rec_voc"))

We can see that the different curves are nearly parallel, which corroborates the lack of significant interaction between the two variables.

Knowing that the interaction between rec_voc and LoE corresponds to the change in the slope of LoE for every one unit increase in rec_voc (or vice-versa), we can assess the significance of this interaction by looking at the contrast / difference between two slopes of rec_voc separated by one unit increase in LoE (the result will be independent from the choice of the two values separated by one unit).

list.LoE.red <- list(LoE = seq(median(df_reduced$LoE) - 0.5, median(df_reduced$LoE) + 0.5, by = 1))
emtrends(model_rep, pairwise ~ LoE, var = "rec_voc", at = list.LoE.red, adjust = "mvt", infer = c(T, T))

$emtrends
 LoE rec_voc.trend     SE  df asymp.LCL asymp.UCL z.ratio p.value
 184        0.0267 0.0203 Inf    -0.013    0.0665   1.320  0.1869
 184        0.0268 0.0203 Inf    -0.013    0.0667   1.320  0.1870

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 

$contrasts
 contrast             estimate      SE  df asymp.LCL asymp.UCL z.ratio p.value
 LoE183.5 - LoE184.5 -0.000101 0.00021 Inf -0.000511   0.00031  -0.481  0.6309

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 

We find that the estimate for the interaction is not significantly different from 0 - the p-value is much larger than 0.05.

For age : phono_mem:

plot_model(model_rep, type = "emm", terms = c("age [all]", "phono_mem"))

list.phono_mem.red <- list(phono_mem = seq(median(df_reduced$phono_mem) - 0.5, median(df_reduced$phono_mem) + 0.5, by = 1))
emtrends(model_rep, pairwise ~ phono_mem, var = "age", at = list.phono_mem.red, adjust = "mvt", infer = c(T, T))

$emtrends
 phono_mem age.trend     SE  df asymp.LCL asymp.UCL z.ratio p.value
       3.5    0.0148 0.0109 Inf  -0.00657    0.0361   1.357  0.1748
       4.5    0.0176 0.0151 Inf  -0.01200    0.0473   1.166  0.2434

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 

$contrasts
 contrast                    estimate     SE  df asymp.LCL asymp.UCL z.ratio p.value
 phono_mem3.5 - phono_mem4.5 -0.00286 0.0125 Inf   -0.0274    0.0217  -0.228  0.8193

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 

Once again, the p-value is much larger than 0.05.

For age : rec_voc:

plot_model(model_rep, type = "emm", terms = c("age [all]", "rec_voc"))

list.rec_voc.red <- list(rec_voc = seq(median(df_reduced$rec_voc) - 0.5, median(df_reduced$rec_voc) + 0.5, by = 1))
emtrends(model_rep, pairwise ~ rec_voc, var = "age", at = list.rec_voc.red, adjust = "mvt", infer = c(T, T))

$emtrends
 rec_voc age.trend     SE  df asymp.LCL asymp.UCL z.ratio p.value
      20    0.0165 0.0111 Inf  -0.00522    0.0382   1.488  0.1367
      21    0.0150 0.0112 Inf  -0.00698    0.0369   1.336  0.1816

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 

$contrasts
 contrast              estimate      SE  df asymp.LCL asymp.UCL z.ratio p.value
 rec_voc20 - rec_voc21  0.00152 0.00142 Inf  -0.00125    0.0043   1.075  0.2826

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 

The p-value for the interaction is higher than 0.05.

For phono_mem : rec_voc:

plot_model(model_rep, type = "emm", terms = c("phono_mem", "rec_voc"))

list.phono_mem.red <- list(phono_mem = seq(median(df_reduced$phono_mem) - 0.5, median(df_reduced$phono_mem) + 0.5, by = 1))
emtrends(model_rep, pairwise ~ phono_mem, var = "rec_voc", at = list.phono_mem.red, adjust = "mvt", infer = c(T, T))

$emtrends
 phono_mem rec_voc.trend     SE  df asymp.LCL asymp.UCL z.ratio p.value
       3.5        0.0207 0.0267 Inf  -0.03161    0.0731   0.776  0.4376
       4.5        0.0370 0.0239 Inf  -0.00988    0.0839   1.547  0.1219

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 

$contrasts
 contrast                    estimate     SE  df asymp.LCL asymp.UCL z.ratio p.value
 phono_mem3.5 - phono_mem4.5  -0.0163 0.0324 Inf   -0.0799    0.0473  -0.502  0.6158

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 

The p-value is once again much larger than 0.05.

Our investigation of the interactions there shows that none of them is statistically significant. A possible option would then be to simplify the model by dropping these interactions. However, this amounts to model selection (for the fixed effects), which is warned against by a number of prominent statisticians. In what follows, we are therefore going to assess main effects despite the presence of interactions in the model.

6.1.2 Investigating the main effects

Given that none of the interactions we thought could be significant appears to be so, we can focus on the main effects in our model, i.e. the effects of the item-related categorical variables occurrence_l, branching_onset, and V, and the subject-related continous variables rec_voc, phono_mem, LoE, age, L1_syll_complexity and phono_awareness.

For occurrence_l:

plot_model(model_rep, type = "emm", terms = "occurrence_l")

summary(emmeans(model_rep, pairwise ~ occurrence_l, adjust = "mvt", side = "<"), infer = c(TRUE, TRUE), null = 0)$contrasts

 contrast      estimate    SE  df asymp.LCL asymp.UCL z.ratio p.value
 coda - final    -0.351 0.400 Inf      -Inf     0.476  -0.877  0.4098
 coda - other    -1.364 0.336 Inf      -Inf    -0.670  -4.061  0.0001
 final - other   -1.013 0.257 Inf      -Inf    -0.483  -3.949  0.0001

Results are averaged over the levels of: branching_onset, V 
Results are given on the log odds ratio (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: mvt method for 3 estimates 
P value adjustment: mvt method for 3 tests 
P values are left-tailed 

We set the parameter side to reflect our hypothesis that a nonword gets easier to repeat when shifting from l in coda position to l in final position to another structure.

We find that:

• nonwords are significantly more difficult to repeat when l appears in internal coda position than when there is no l (coda - other, p < 0.0001).
• nonwords are significantly more difficult to repeat when l appears in final position than when there is no l (final - other, p < 0.0001)
• nonwords are not significantly more difficult to repeat when l appears in internal coda position than when l appears in final position (coda - final, p = 0.4101)

For branching_onset:

plot_model(model_rep, type = "emm", terms = "branching_onset")

summary(emmeans(model_rep, pairwise ~ branching_onset, adjust = "mvt", side = ">"), infer = c(TRUE, TRUE), null = 0)$contrasts

 contrast                            estimate    SE  df asymp.LCL asymp.UCL z.ratio p.value
 branching_onset0 - branching_onset1    0.812 0.172 Inf    0.4654       Inf   4.725  <.0001
 branching_onset0 - branching_onset2    1.830 0.465 Inf    0.8928       Inf   3.939  0.0001
 branching_onset1 - branching_onset2    1.018 0.471 Inf    0.0678       Inf   2.162  0.0358

Results are averaged over the levels of: occurrence_l, V 
Results are given on the log odds ratio (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: mvt method for 3 estimates 
P value adjustment: mvt method for 3 tests 
P values are right-tailed 

We set the parameter side to reflect our hypothesis that the more branching onsets in a nonword, the more difficult it is to repeat.

We observe that all the contrasts are significant, and that therefore:

• A nonword is more difficult to repeat when it has 1 branching onset than when it has none (p < 0.0001)
• A nonword is more difficult to repeat when it has 2 branching onsets than when it has none (p < 0.0001)
• A nonword is more difficult to repeat when it has 2 branching onsets than when it has 1 (p = 0.0355)

For V:

plot_model(model_rep, type = "emm", terms = "V")

summary(emmeans(model_rep, pairwise ~ V, adjust = "mvt", side = ">"), infer = c(TRUE, TRUE), null = 0)$contrasts

 contrast estimate    SE  df asymp.LCL asymp.UCL z.ratio p.value
 V1 - V2     0.101 0.199 Inf   -0.3133       Inf   0.508  0.5926
 V1 - V3     0.595 0.203 Inf    0.1735       Inf   2.939  0.0048
 V2 - V3     0.494 0.201 Inf    0.0765       Inf   2.464  0.0190

Results are averaged over the levels of: occurrence_l, branching_onset 
Results are given on the log odds ratio (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: mvt method for 3 estimates 
P value adjustment: mvt method for 3 tests 
P values are right-tailed 

The parameter side corresponds to the hypothesis that the less vowels in a nonword, the easier it is to repeat.

We find two significant differences: a nonword with a single vowel is easier to repeat than a nonword with 3 vowels (p = 0.0047), and a nonword with 2 vowels is easier to repeat than a nonword with 3 vowels (p = 0.0192)

For rec_voc:

plot_model(model_rep, type = "emm", terms = "rec_voc [all]")

summary(emtrends(model_rep, ~ rec_voc, var = "rec_voc", adjust = "mvt", side = ">"), infer = c(TRUE, TRUE), null = 0)

 rec_voc rec_voc.trend     SE  df asymp.LCL asymp.UCL z.ratio p.value
    20.4        0.0268 0.0203 Inf  -0.00662       Inf   1.320  0.0935

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 
P values are right-tailed 

We set the parameter side to reflect our hypothesis that the larger a child’s French receptive vocabulary, higher the probability of correct repetition.

We observe only a tendency for rec_voc (p = 0.0935), suggesting that the size of a child’s receptive vocabulary positively impacts her/his ability to correctly repeat nonwords.

For phono_mem:

plot_model(model_rep, type = "emm", terms = "phono_mem")

summary(emtrends(model_rep, ~ phono_mem, var = "phono_mem", adjust = "mvt", side = ">"), infer = c(TRUE, TRUE), null = 0)

 phono_mem phono_mem.trend    SE  df asymp.LCL asymp.UCL z.ratio p.value
      3.88          0.0521 0.126 Inf    -0.155       Inf   0.415  0.3392

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 
P values are right-tailed 

We set the parameter side to reflect our hypothesis that the larger a child’s phonological memory, the higher the probability of correct repetition.

We do not observe a significant effect of phono_mem.

For LoE:

plot_model(model_rep, type = "emm", terms = "LoE [all]")

summary(emtrends(model_rep, ~ LoE, var = "LoE", adjust = "mvt", side = ">"), infer = c(TRUE, TRUE), null = 0)

 LoE LoE.trend      SE  df asymp.LCL asymp.UCL z.ratio p.value
 184   0.00182 0.00159 Inf -0.000786       Inf   1.149  0.1253

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 
P values are right-tailed 

We set the parameter side to reflect our hypothesis that the longer the exposure to French, the higher the probability of correct repetition.

We do not observe a significant effect of LoE on the probability of correct repetition.

For age:

plot_model(model_rep, type = "emm", terms = "age [all]")

summary(emtrends(model_rep, ~ age, var = "age", adjust = "mvt", side = ">"), infer = c(TRUE, TRUE), null = 0)

  age age.trend     SE  df asymp.LCL asymp.UCL z.ratio p.value
 90.5    0.0158 0.0111 Inf  -0.00241       Inf   1.428  0.0766

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 
P values are right-tailed 

We set the parameter side to reflect our hypothesis that the older a child is, the higher the probability of correct repetition.

We do not observe a significant effect of age, although we find a tendency suggesting that the older a child, the higher the probability of correct repetition (p = 0.0766).

For L1_syll_complexity:

plot_model(model_rep, type = "emm", terms = "L1_syll_complexity [all]")

summary(emtrends(model_rep, ~ age, var = "L1_syll_complexity", adjust = "mvt", side = ">"), infer = c(TRUE, TRUE), null = 0)

  age L1_syll_complexity.trend    SE  df asymp.LCL asymp.UCL z.ratio p.value
 90.5                   0.0389 0.115 Inf    -0.151       Inf   0.337  0.3681

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 
P values are right-tailed 

We set the parameter side to reflect our hypothesis that the higher the complexity of the syllables of the L1, the higher the probability of correct repetition.

We do not observe a significant effect of L1_syll_complexity (p = 0.3681).

For phono_awareness:

plot_model(model_rep, type = "emm", terms = "phono_awareness [all]")

summary(emtrends(model_rep, ~ phono_awareness, var = "phono_awareness", adjust = "mvt", side = ">"), infer = c(TRUE, TRUE), null = 0)

 phono_awareness phono_awareness.trend    SE  df asymp.LCL asymp.UCL z.ratio p.value
           -2.82                 0.161 0.068 Inf    0.0494       Inf   2.371  0.0089

Results are averaged over the levels of: occurrence_l, branching_onset, V 
Confidence level used: 0.95 
P values are right-tailed 

We set the parameter side to reflect our hypothesis that the more developed a child’s phonological awareness, the higher the probability of correct repetition.

We observe a significant effect of phono_awareness (p = 0.0089): the more developed the child’s phonological awareness, the better they are at correctly repeating the nonwords.

6.1.3 Summary for the main effects

All our hypotheses were not confirmed. We found the following results:

• nonwords are significantly more difficult to repeat when l appears in internal coda position than when there is no l (p < 0.0001).
• nonwords are significantly more difficult to repeat when l appears in final position than when there is no l (p < 0.0001).
• the more branching onsets nonwords have, the more difficult they are to repeat (1 versus 0 BO: p < 0.0001; 2 versus 0 BO: p < 0.0001; 2 versus 1 BO: p = 0.0355)
• nonwords with 3 vowels are more difficult to repeat than nonwords with 1 or 2 vowels (p = 0.0047 and p = 0.0192, respectively)
• the more developed a child’s phonological awareness, the easier for them to repeat the nonwords (p < 0.0089)

Additionally, we observed two statistical tendencies:

• the larger a child’s French receptive vocabulary, the easier for them to repeat the nonwords (p < 0.0935)
• the older a child, the easier for them to repeat the nonwords (p = 0.0766)

We did not observe effects of phono_mem or of LoE. We also did not find a significant difference between l appearing in internal coda position in nonwords and l appearing in final position. We finally did not find an effect of the complexity of syllables in L1.

Our analysis overall proves to be quite simple, in the sense that we did not observe any significant interaction.

p1 <- plot_model(model_rep, type = "emm", terms = "occurrence_l", show.values = TRUE, value.offset = .3, colors = "gs") +
  labs(title = "Effect of occurence_l", x = "Location of /l/", y = "Predicted probability of correct repetition")
p2 <- plot_model(model_rep, type = "emm", terms = "branching_onset", show.values = TRUE, value.offset = .3, colors = "gs") +
  labs(title = "Effect of branching_onset", x = "Number of branching onsets", y = "Predicted probability of correct repetition")
p3 <- plot_model(model_rep, type = "emm", terms = "V", show.values = TRUE, value.offset = .3, colors = "gs") +
  labs(title = "Effect of V", x = "Number of vowels", y = "Predicted probability of correct repetition")
p4 <- plot_model(model_rep, type = "emm", terms = "rec_voc [all]", show.values = TRUE, value.offset = .3, colors = "gs") +
  labs(title = "Effect of rec_voc", x = "Size of the receptive vocabulary", y = "Predicted probability of correct repetition")
p5 <- plot_model(model_rep, type = "emm", terms = "age [all]", show.values = TRUE, value.offset = .3, colors = "gs") +
  labs(title = "Effect of age", x = "Age", y = "Predicted probability of correct repetition")
p6 <- plot_model(model_rep, type = "emm", terms = "phono_awareness [all]", show.values = TRUE, value.offset = .3, colors = "gs") +
  labs(title = "Effect of phonological awareness", x = "Phonological awareness", y = "Predicted probability of correct repetition")

grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 3)

6.2 Analyzing the random effects

While our analysis is centered on fixed effects, it is interesting to pay attention to the distribution of values of the different random effects in our model.

p <- plot_model(model_rep, type="re", sort.est="sort.all", grid=F, transform = NULL)

The most interesting source of information is likely the distribution of the random intercepts of the different children’s L1:

p[[3]] + scale_color_grey()

We can notice that the languages closests to French genealogically speaking, i.e., Spanish, Italian, Portuguese and Brazilian Portuguese have negative intercepts (while accounting for the complexity of their syllables with L1_syll_complexity as a fixed continuous predictor). Romanian, however, has the second highest positive intercept. The distribution does not seem suggestive of an effect of earlier bilingualism (before the acquisition of French) - the intercepts for Arabic/Italian, Russian/Chechen etc. are rather in the middle of the distribution.

p[[1]] + scale_color_grey()

There is no striking element in the distribution observed.

Finally, we can display the values of the random intercepts for subjects:

p[[2]] + scale_color_grey()

Once again, there is no striking element in the distribution observed.